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ABSTRACT 

The  dynamical  evolution  of  terrestrial  planets  resembling  Mercury  in  the  vicinity  of  spin-orbit  resonances  is 
investigated  using  comprehensive  harmonic  expansions  of  the  tidal  torque  taking  into  account  the  frequency- 
dependent  quality  factors  and  Love  numbers.  The  torque  equations  are  integrated  numerically  with  a  small  step 
in  time,  including  the  oscillating  triaxial  torque  components  but  neglecting  the  layered  structure  of  the  planet  and 
assuming  a  zero  obliquity.  We  find  that  a  Mercury-like  planet  with  a  current  value  of  orbital  eccentricity  (0.2056) 
is  always  captured  in  3:2  resonance.  The  probability  of  capture  in  the  higher  2:1  resonance  is  approximately  0.23. 
These  results  are  confirmed  by  a  semi-analytical  estimation  of  capture  probabilities  as  functions  of  eccentricity  for 
both  prograde  and  retrograde  evolutions  of  spin  rate.  As  follows  from  analysis  of  equilibrium  torques,  entrapment 
in  3:2  resonance  is  inevitable  at  eccentricities  between  0.2  and  0.41.  Considering  the  phase  space  parameters  at  the 
times  of  periastron,  the  range  of  spin  rates  and  phase  angles  for  which  an  immediate  resonance  passage  is  triggered 
is  very  narrow,  and  yet  a  planet  like  Mercury  rarely  fails  to  align  itself  into  this  state  of  unstable  equilibrium  before 
it  traverses  2:1  resonance. 
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1.  INTRODUCTION 

Planets  of  terrestrial  type  with  solid  mantles  are  subject  to 
triaxial  and  tidal  torques  exerted  by  their  host  star.  A  planet 
just  formed  from  the  protoplanetary  disk  rotates  at  a  much 
higher  rate  than  the  rate  of  its  orbital  motion.  In  the  course 
of  millions  to  billions  of  years,  the  secular  terms  of  the  tidal 
torque  cause  the  planet  to  spin  down.  The  tidal  bulges  moving 
across  the  planet  at  different  frequencies  result  in  a  gradual 
loss  of  kinetic  energy  through  friction  and  heating.  The  energy 
dissipation  rate  is  normally  so  slow  that  most  of  the  major  planets 
in  the  solar  system  still  rotate  faster  than  they  revolve  around 
the  Sun,  with  the  exception  of  Venus  with  its  slow  retrograde 
rotation  and  Mercury,  which  is  in  a  3:2  spin-orbit  resonance 
(Pettengill  &  Dyce  1965).  Presumably,  the  planet  traversed  a 
number  of  higher-order  resonances  before  it  reached  this  state. 
The  ultimate,  and  the  most  stable,  state  for  a  rotating  planet 
subject  to  tidal  forces  is  1 : 1  resonance,  when  the  rotation  rate  is 
equal  to  the  orbital  rate,  and  the  planet  is  always  pointing  with 
its  most  elongated  dimension  toward  the  star.  The  dissipation 
of  the  energy  of  rotation  also  diminishes  the  obliquity  of  the 
planet’s  equator,  gradually  aligning  the  axes  of  rotation  and  of 
the  orbit.  Mercury’s  dynamical  evolution  has  been  much  faster 
than  that  of  other  solar  planets,  because  (1)  it  is  closer  to  the  Sun 
and  (2)  its  orbital  eccentricity  varied  in  a  relatively  wide  range 
and  has  been  higher  than  the  eccentricity  of  the  other  close-in 
planets.  The  importance  of  eccentricity  for  Mercury’s  chaotic 
evolution  was  emphasized  by  Correia  &  Laskar  (2004,  2009), 
who  upwardly  revised  the  original  estimate  of  the  probability  of 
capture  to  the  present  3:2  resonance  at  7%  (Goldreich  &  Peale 
1966). 

Mercury  can  serve  as  a  good  model  for  smaller-mass,  rocky 
exoplanets,  especially  those  orbiting  M  dwarfs  in  their  habitable 
zones.  Although  the  exact  rheology  of  exoplanets  will  remain 
a  matter  of  speculation  for  the  foreseeable  future,  it  seems 
reasonable  to  adopt  the  parameters  and  models  obtained  for  the 


Earth  and  the  Moon.  The  objective  of  this  paper  is  to  investigate 
the  circumstances  of  the  transition  of  a  Mercury-like  planet 
with  an  Earth-like  rheology  through  a  high-order  spin-orbit 
resonance.  In  this  case,  the  widely  accepted  approximations 
for  the  value  of  tidal  torque  in  the  vicinity  of  a  resonance 
are  not  applicable.  Eurthermore,  the  oscillatory  terms  of  the 
force  cannot  be  neglected.  We  employ  high-order  expansions 
of  the  torque  in  harmonics  of  tidal  frequency  and  powers  of 
eccentricity,  and  a  relation  for  the  Love  number  as  a  function 
of  tidal  frequency  in  terms  of  real  and  imaginary  compliances 
(Section  2).  The  resulting  differential  equation  of  second  order, 
which  includes  both  the  tidal  and  triaxial  torque  components,  is 
integrated  with  a  step  much  smaller  than  the  period  of  rotation, 
with  the  current  best  estimates  for  Mercury  (Section  3).  With 
the  current  value  of  eccentricity  (e  =  0.20563),  the  planet 
traverses  2:1  resonance  with  an  estimated  probability  of  0.77, 
and  is  always  captured  in  3:2  resonance.  The  transition  is  very 
fast,  accompanied  by  a  significant  step-down  in  the  average 
rotation  rate.  In  Section  4,  more  integrations  are  performed  with 
the  initial  rate  set  to  exactly  a  2: 1  resonance  and  with  various 
initial  phase  angles.  It  is  revealed  that  such  a  planet  inserted  in 
resonance  almost  always  stays  in  it.  The  actual  passage  through 
resonance  can  only  occur  through  a  tiny  area  of  the  phase  space 
around  a  phase  angle  of  H-7r/2  or  — 7r/2.  A  qualitative  physical 
explanation  for  this  curious  fact  is  given.  We  discuss  in  the 
final  Section  5  why  our  model  planet  is  unlikely  to  traverse 
2:1  resonance  despite  the  tightness  of  the  required  conditions 
to  do  so.  We  derive  the  probabilities  of  capture  in  spin-orbit 
resonances  for  both  prograde  and  retrograde  evolutions  of 
spin  rate,  and  compare  them  with  the  ranges  of  equilibrium 
eccentricities. 

2.  HARMONIC  EXPANSIONS  OE  TORQUES 

The  instantaneous  torque  acting  on  a  rotating  planet  is 
the  sum  of  the  triaxial  torque,  caused  by  the  quadrupole 
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Table  1 

Explanation  of  Notations 


Notation 

Description 

? 

Moment  of  inertia  coefficient 

R 

Radius  of  planet 

T 

Torque 

M2 

Mass  of  planet 

Ml 

Mass  of  star 

a 

Semimajor  axis  of  planet 

r 

Instantaneous  distance  of  planet  from  star 

V 

True  anomaly  of  planet 

e 

Orbital  eccentricity 

M 

Mean  anomaly  of  planet 

B 

Second  moment  of  inertia 

A 

Third  moment  of  inertia 

C 

Moment  of  inertia  around  spin  axis 

n 

Mean  motion,  i.e.,  27r/Porb 

Q 

Gravitational  constant,  =  66468  m^  kg~^  yr“^ 

t:m 

Maxwell  time,  i.e.,  ratio  of  viscosity  to  unrelaxed  rigidity^ 

Unrelaxed  rigidity  modulus 

a 

Tidal  lag  responsivity 

Note.  ^  For  the  differences  between  relaxed  and  unrelaxed  rigidity  moduli,  see 
Efroimsky  (2012a,  2012b). 


inertial  momentum,  and  the  tidal  torque,  caused  by  the  dynamic 
deformation  of  its  body.  In  neglect  of  the  obliquity  (Danby 
1962), 


TtrI  +  TflDE 


(1) 


with  9  being  the  sidereal  angle  of  the  planet  reckoned  from  the 
axis  of  its  largest  elongation.  All  other  notations  used  in  this 
formula  and  throughout  the  paper  are  explained  in  Table  1 .  We 
consider  the  specific  but  representative  case  when  the  obliquity 
of  the  planet’s  equator  is  small  (i  ~  0)  and  the  planet  is  not 
too  close  to  the  star  {R/a  1).  Neglecting  the  precession 
and  nutation  of  the  planet,  the  triaxial  torque  is  (Danby  1962; 
Goldreich  &  Peale  1966) 


Ttri  — 


3 

— {B  —  A)n^—  sin  2(0  —  v). 
2 


(2) 


Using  the  comprehensive  development  of  Kaula  and  Darwin’s 
harmonic  decomposition  of  the  tidal  torque  by  Efroimsky  & 
Lainey  (2007),  Efroimsky  &  Williams  (2009),  and  Efroimsky 
(2012a,  2012b),  one  can  write  a  simplified  equation  for  the  tidal 
torque,  which  we  call  Efroimsky’s  torque 


3  R^  ^  ^ 

Ttide  =  ^  G20qie)  ^  G20j(e) 

^  9=-i  i=-i 

X  [Kc(2,  X220q)Signico220q)cos{(q  -  j)M) 

+  Ks{2,  X220q)  sin  {{q  -  ;)M)],  (3) 

where  the  positively  defined  forcing  frequency 

X22Qq  =  l®220^l  =  l(2  +  ^)n  —  20|,  (4) 

the  so-called  G-functions  of  Kaula  G20 j  are  related  to  power  se¬ 
ries  in  eccentricity  via  Hansen’s  coefficients  (e.g.,  Dobrovolskis 
1995): 

G2Qj{e)  ^  (5) 


and  the  all-important  “quality  functions”  are 


KcU,  X)  = 


3 _ A/x^(x) 

2(/-  l)(lR(x)  +  A,x)^+A2(x) 


(6) 


^  3  lii^(x)  +  a"(x)  +  Axiii(x) 

2(/-l)  (Jli(x)+A,x)2+32(x) 

Einally,  the  remaining  terms  in  this  equation  are 

4;r(2/2  H- 4/ -h  3)/?V 

A;  =  - ~ - , 

3igMl 

and 


(7) 


(8) 


S1'i(x)  =  X+x'  "r^“cos(^^r(l  H-a) 

A(X)  =  -  -  X  sin  (^)  r(l  +  a).  (9) 

The  lli(x)  and  Ss(x)  functions  are  the  real  and  imaginary  parts 
of  the  complex  compliance,  respectively.  Equation  (3)  includes 
oscillating  terms  of  the  tidal  torque,  q  ^  j.  Eor  a  planet  similar 
to  Mercury,  these  oscillating  terms  prove  very  small  and  can  be 
safely  neglected.^  The  results  of  integration  over  20,000  years 
were  hardly  different  with  or  without  these  terms.  Omitting  the 
oscillating  terms  of  the  tidal  torque,  we  arrive  at  the  simplified 
expression  for  the  secular  part  of  the  torque 

3  R^ 

Ttide  =  ^2  ^20q(^')^c(2,  X220q)Sign{a>220q)- 

^  ?=-i 

(10) 

Note  that  expansions  (3)  and  (10)  are  limited  to  /  =  2. 
The  higher-order  expansion  terms  (/  =  3)  are  neglected  be¬ 
cause  the  coefficients  are  at  least  eight  orders  of  magnitude 
smaller  than  the  1  —  2  terms  for  Mercury  or  similar  plan¬ 
ets.  The  expansion  in  harmonics  of  mean  anomaly  is  limited 
to  the  range  q,  j  —  —  1 ,  0,  . . . ,  4  for  the  ease  of  compu¬ 
tation,  the  other  terms  being  much  smaller  than  these  six. 
Eor  example,  G205(0. 20563)  =  0.019,  to  be  compared  with 
G2oo(0. 20563)  =  0.896,  and  the  rest  of  the  Kaula  coefficients 
for  the  omitted  terms  are  even  smaller.  The  crucial  differences 
between  this  approach  and  the  previous  studies,  e.g.,  by  Celletti 
et  al.  (2007),  is  that  more  realistic  relations  are  employed  here 
for  the  “quality  functions”  of  tidal  frequency  (Equations  (6)  and 
(7))  with  rheological  parameters  supported  by  observations  and 
theory,  and  that  full  account  is  taken  of  the  fast,  oscillating  terms 
of  the  triaxial  torque,  which  are  integrated  in  a  small  step. 

3.  INTEGRATION  OE  MERCURY 

The  evolution  of  Mercury’s  spin  can  be  considered  as  a 
template  model  of  a  rocky,  relatively  small  planet  whose 
eccentricity  can  be  driven  to  sufficiently  large  values  due  to  the 
interactions  with  larger  planets  in  the  system.  It  is  also  one  of 
the  few  planets  whose  physical  parameters  are  fairly  well  known 
(Anderson  et  al.  1987).  The  parameters  used  in  our  computations 
are  listed  in  Table  2.  The  Maxwell  time  Xm  was  set  to  Earth’s 
value,  assuming  a  similar  composition  and  temperature.  The 
sensitivity  of  this  crucial  parameter  to  temperature  has  to  be 

*  The  oscillating  components  of  the  tidal  torque  cannot  be  neglected  for  a 
nearly  spherical  planet,  i.e.,  when  (B  —  A)/C  is  very  small. 
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Figure  1.  Rotation  acceleration  caused  by  the  secular  tidal  torque  (i.e.,  j  =  q  in  Equation  (3),  or  Equation  (10))  on  a  Mercury-like  planet  in  the  vicinity  of  2:1 
resonance  (left)  and  3:2  resonance  (right).  Dramatic  variations  of  the  torque  are  confined  to  very  narrow  ranges  of  frequencies,  while  the  effective  torque  integrated 
over  a  wider  range  of  frequencies  is  negative.  The  variable,  antisymmetric  part  of  the  torque  is  caused  by  the  given  resonance.  The  almost  constant  negative  offset  is 
due  to  the  contribution  of  other  resonances.  Despite  their  appearance  on  the  chosen  scale,  the  functions  are  smooth  and  differentiable. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


Table  2 

Parameters  of  the  Mercury  Integration 


Parameter 

Value 

? 

2 

5 

R 

2.44  X  10®  m 

M2 

3.3  X  10^3  kg 

Ml 

1.99  X  103°  kg 

a 

5.79  X  10‘“  m 

e 

0.20563 

{B  -  A)/C 

1.2  X  lO'"* 

n 

26.088  yr-' 

rM 

500  yr 

0.8  X  10"  kgm-'  s-3 

a 

0.2 

taken  into  account  for  rocky  exoplanets  whose  interiors  may 
be  significantly  hotter.  Higher  temperatures  result  in  smaller 
viscosity  and,  thus,  much  shorter  Maxwell  times.  The  rigidity 
/r  was  also  set  to  the  value  for  the  Earth.  The  characteristic  time 
ta  is  in  fact  frequency-dependent,  and  we  modeled  it  as 

taCx)  =  50000  exp(-x/0.2)-H  500  (H) 

in  years.  This  expression  is  chosen  to  represent  the  expected 
behavior  of  r^,  which  sharply  rises  toward  small  frequencies. 
However,  the  results  proved  quite  insensitive  to  the  functional 
form  of  ta-  In  fact,  practically  the  same  conclusions  can  be 
obtained  by  simply  setting  ta  =  Tm-  Again,  this  may  not  be  the 
case  for  hotter,  less  viscous  planets. 

The  tidal  torque  at  a  given  eccentricity  is  a  slowly  varying 
function  of  spin  rate  everywhere  except  the  vicinity  of  spin- 
orbit  resonances  0  =  (1  h-  ql2)n.  Figure  1  shows  in  detail  the 
dependence  of  the  overall  angular  acceleration  6  of  the  planet 
caused  by  the  secular  tidal  torque  only,  in  the  vicinity  of  2:1 
{q  =  2)  and  3:2  (q  =  1)  resonances.  Note  that  the  dramatic 
changes  of  tidal  torque  take  place  within  a  very  narrow  interval 
of  tidal  frequencies.  The  widely  used  assumption  that  0  oc  —6 
is  justified  only  in  a  vanishingly  small  range  of  spin  rates,  i.e., 
6  e  [1.9999977,  2.0000023]  n  in  the  case  of  2:1  resonance  and 
for  the  chosen  value  of  the  Maxwell  time  Tm-  The  amplitude 
of  the  oscillatory  triaxial  torques  outside  the  1:1  resonance  is 
a  few  orders  of  magnitude  greater  than  the  peak  values  of  the 
tidal  torque  (Figure  2).  The  scale  of  free  librations,  caused  by 
the  triaxial  torque,  is  also  a  few  orders  of  magnitude  greater  than 


0.0  0.5  1.0  1.5  2.0 


time  years 

Figure  2.  Rotation  acceleration  caused  by  the  triaxial  torque  on  a  Mercury-like 
planet  as  a  function  of  time  at  ^  =  1.6  n. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 

the  characteristic  width  of  the  kink  in  Figure  1 .  The  net  effect 
of  the  secular  tidal  force  integrated  over  one  libration  period 
is  to  spin  the  planet  down.  This  follows  from  the  fact  that  the 
mean  value  of  acceleration  at  0  =  2n  integrated  over  a  much 
wider  interval  than  the  width  of  the  kink,  is  negative  for  these 
resonances.  The  peak  torque  below  the  resonance  rate  tapers  off 
quickly  and  becomes  negative  at  0  =  1.999895  n.  The  torque 
function  near  3:2  resonance  (Figure  1,  right)  is  similar  in  shape, 
but  the  amplitude  is  larger  and  the  positive  shoulder  below  the 
resonance  is  much  broader.  The  1 : 1  resonance  is  fundamentally 
different  from  the  higher-order  counterparts  in  that  the  tidal 
torque  is  positive  at  all  rotation  rates  between  0  and  n. 

3.1.  Probabilities  of  Capture  in  2:1  and  3:2  Resonances 

The  ordinary  differential  equation  (1)  was  integrated  with  a 
grid  of  initial  conditions  0(0),  0(0)  for  30,000-55,000  years 
with  a  maximum  step  of  2  x  10“^  yr.  The  initial  mean  anomaly 
was  always  set  to  M(0)  =  0.  For  2:1  resonance,  a  grid  of 
40  integrations  with  0(0)  =  2.013 n  and  0(0)  =  7r//40, 
i  =  0,  1, . . . ,  39  was  performed.  These  integrations  resulted 
in  9  captures  and  3 1  passages  of  the  resonance.  Therefore,  the 
estimated  probability  of  2: 1  capture  is  0.23  at  the  current  value  of 
eccentricity.  The  amplitude  of  free  libration  gradually  increases 
as  the  planet  spins  down  toward  the  point  of  resonance.  Once  a 
lower  swing  in  0  reaches  the  resonance,  a  fast  transition  into  a 
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new  spin  state  occurs  on  a  timescale  of  a  few  years.  The  average 
rate  takes  a  sudden  leap  to  a  significantly  smaller  value,  and  the 
libration  starts  to  decrease.  This  set  of  simulations  shows  that 
the  chance  of  Mercury  being  captured  into  the  2: 1  resonance  at 
the  current  value  of  eccentricity  is  modest. 

The  outcome  of  integrations  in  the  vicinity  of  3:2  resonance  is 
quite  different.  A  similar  set  of  integrations  with  0(0)  =  1.518n 
and0(O)  =  it  i/AQ,i  =  0,  1,  . . . ,  39  resulted  in 40 captures.  The 
estimated  probability  of  3:2  capture  is  1.  As  soon  as  9  reaches 
the  point  of  resonance,  the  amplitude  of  libration  abruptly 
doubles,  but  the  mean  rotation  stays  around  1.5  n.  Once  the 
planet  is  captured,  the  amplitude  of  libration  starts  to  slowly 
decline.  The  period  of  libration  starts  to  decrease  immediately 
after  the  capture.  A  longer  integration  for  100,000  yr  shows 
that  the  period  of  libration  asymptotically  approaches  16  yr, 
which  is  the  theoretical  value  (Murray  &  Dermott  2000)  for  the 
physical  parameters  listed  in  Table  2.  This  confirms  the  validity 
of  our  computations  of  the  triaxial  torque.  Further  numerical 
experiments  revealed  that  the  planet  traverses  3:2  resonance  in 
a  similar  manner  to  2: 1  resonance  at  significantly  smaller  values 
of  eccentricity.  We  conclude  that  our  numerical  simulations 
are  consistent  with  the  fact  that  Mercury  is  entrapped  in  a 
3:2  resonance  with  the  current  value  of  eccentricity.  This 
outcome  is  indeed  most  likely  unless  the  eccentricity  acquired 
much  different  values  in  the  past.  An  even  higher  eccentricity 
would  have  resulted  in  entrapment  in  a  2:1  resonance  long  in 
the  past  when  the  planet  was  approaching  this  rotation  rate. 
Alternatively,  a  considerably  smaller  eccentricity  at  the  moment 
of  approaching  3:2  resonance  would  have  made  the  planet 
traverse  it  quickly  and  continue  to  spin  down. 

3.2.  Equilibrium  Torques  and  Evolution  of  Eccentricity 

Correia  &  Laskar  (2004)  pointed  out  that  Mercury’s  eccen¬ 
tricity  has  varied  chaotically  during  its  long  dynamical  evolution 
in  the  solar  system.  The  current  state  of  Mercury,  entrapped  in 
a  3:2  spin-orbit  resonance,  is  the  result  of  a  long  history  of  tidal 
and  orbital  interactions,  probably  marked  by  multiple  passages 
of  spin-orbit  resonances.  In  particular,  if  the  initial  spin  rate 
of  Mercury  was  much  greater  than  it  is  today,  the  planet  has 
successfully  traversed  a  number  of  higher  resonances.  Under¬ 
standing  the  dependence  of  tidal  torque  on  eccentricity  in  the 
framework  of  Efroimsky’s  model  is  as  important  as  the  truthful 
estimation  of  capture  probabilities  in  each  resonance. 

Correia  &  Laskar  (2004,  2009)  employed  a  linear  torque 
model,  also  known  as  MacDonald’s  torques  (Goldreich  & 
Peale  1966)  or  the  Constant  Phase  Lag  model  (Ferraz-Mello 
et  al.  2008),  in  which  the  secular  torque  is  linearly  dependent 
on  the  time  derivative  of  tidal  frequency.  As  we  noted  in 
Section  3,  Efroimsky’s  torque  is  linear  only  in  vanishingly 
narrow  intervals  around  resonances  (Ligures  1(a)  and  (b)). 
The  widely  used  assumption  that  the  secular  tidal  torque 
is  linear  everywhere  across  the  range  of  relative  spin  rates 
O/n  leads  to  well-known  problems  and  inconsistencies.  Lor 
example,  a  slowing  down  Moon  in  the  absence  of  quadrupole 
momentum  is  not  allowed  to  descend  into  1:1  resonance 
(synchronous  rotation)  for  any  nonzero  eccentricity  (Murray  & 
Dermott  2000;  Williams  &  Efroimsky  2012),  because  the  linear 
torque,  monotonously  increasing  with  growing  eccentricity, 
changes  sign  from  negative  to  positive  at  a  certain  equilibrium 
eccentricity  for  any  fixed  9/n  >  1.  Therefore,  a  perfectly 
spherical  Moon  should  have  stalled  in  its  de-spinning  at  a  higher 
than  synchronous  rate,  namely,  0equ/«  =  \+6e^  =  1.018. 
This  is  not  the  case  for  Efroimsky’s  tidal  torque,  which  sharply 


Figure  3.  Equilibrium  eccentricity  of  a  Mercury-like  planet  separating  the  areas 
of  negative  (spin-down)  and  positive  (spin-up)  secular  tidal  torque  for  a  grid  of 
relative  rates  of  rotation.  The  smooth  curve  represents  the  linear  torque  widely 
used  in  the  literature.  The  jagged  curve  represents  Efroimsky’s  torque. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 

decreases  at  supersynchronous  rotation  rates,  counteracting  the 
effects  of  eccentric  motion.  The  dependence  of  equilibrium 
eccentricity,  whence  the  linear  torque  disappears,  for  a  range 
of  rotation  rates  is  depicted  in  Eigure  3  with  the  monotonously 
rising  dotted  line.  It  implies  that  synchronous  rotation,  so 
commonly  observed  among  the  satellites  of  the  solar  system, 
is  not  attainable  for  any  nonzero  eccentricity.  If  Mercury 
traversed  3:2  resonance  and  continued  to  spin  down  toward 
synchronization,  it  would  be  stuck  at  an  equilibrium  rate  of 
1 .24  n  with  the  current  value  of  eccentricity. 

The  character  of  equilibrium  torques  is  profoundly  different 
with  Efroimsky’s  model,  depicted  with  the  jagged  dotted  curve 
in  Eigure  3.  The  curve  was  obtained  by  finding  the  roots  of 
Equation  (10)  in  e  for  a  grid  of  6 /n.  The  grid  resolution  is 
chosen  such  that  only  the  slow  variations  of  the  equilibrium 
eccentricity  outside  the  resonances  are  faithfully  represented. 
The  first  jump  from  Cequ  =  Oat0/«  =  1  to  Cequ  ^  0.29 
just  slightly  above  the  synchronous  rate  implies  that  Mercury 
is  allowed  to  be  captured  in  1:1  resonance  if  the  eccentricity 
is  not  too  large  (e  <  0.29),  evolving  either  from  slower  or 
faster  rotation  rates.  Likewise,  entrapment  in  3:2  resonance  is 
inevitable  for  0.2  <  e  <  0.41  if  Mercury  reaches  this  rate  of 
rotation  from  either  direction.  This  is  consistent  with  the  results 
of  numerical  simulations  in  Section  3.1  that  the  probability  of 
3:2  capture  with  the  current  value  of  eccentricity  is  1.  What 
happens  if  Mercury  reaches  O/n  —  1.5  with  an  eccentricity 
below  0.2?  Capture  is  still  possible,  but  the  probability  declines 
with  decreasing  eccentricity  (see  Section  5).  Entrapment  of  the 
planet  in  2:1  resonance  is  inevitable  for  0.34  <  e  <  0.49,  etc. 
The  teeth  of  the  equilibrium  torque  curve  act  as  very  efficient 
resonance  traps  for  a  planet  like  Mercury,  whose  eccentricity 
varies  in  a  fairly  wide  range  over  billions  of  years  of  dynamical 
evolution. 

4.  CONDITIONS  OF  TRAVERSING  A  RESONANCE 

The  purpose  of  our  next  set  of  simulations  is  to  find  out  under 
which  circumstances  the  planet  traverses  a  spin-orbit  resonance. 
We  have  established  that  the  test  planet  passes  the  point  6  —  2n 
quite  quickly,  within  several  years.  At  some  moment  during  this 
passage,  the  planet  is  at  periastron  with  a  rotation  rate  close  to 
2n.  What  is  the  rotation  angle  9  at  this  time,  and  does  this  value 
matter  for  the  way  this  dynamical  transformation  unfolds?  One 
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Figure  4.  Results  of  integration  of  the  Mercury-like  test  planet  with  initial 
conditions  of  resonance,  ^(0)  =  2n,  0(0)  =  0.  The  spin  rate  remains  in 
resonance  with  the  orbital  motion  indefinitely  long. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


Figure  5.  Similar  to  Figure  4,  but  with  a  different  initial  phase  angle,  0(0)  =  :7r /2. 
The  planet  traverses  the  resonance  after  a  couple  of  libration  oscillations,  in  less 
than  100  yr. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


hundred  integrations  were  conducted  with  the  planet  initially 
already  at  resonance  spin  rate  (0(0)  =  2n),  but  with  initial 
phase  angles  ranging  from  — 7r/2  to  +jt/2  in  equal  steps.  For 
most  of  the  interval  of  possible  initial  angles,  the  planet  is  clearly 
captured  in  resonance,  i.e.,  the  spin  rate  continues  to  oscillate 
around  the  resonant  value.  A  typical  example  of  such  entrapment 
with  the  initial  conditions  0(0)  —  2n,  0(0)  =  0,  is  shown  in 
Figure  4.  The  rate  of  rotation  oscillates  around  the  resonant  point 
roughly  between  1.9996  n  and  2.0006 «,  seemingly  forever. 
The  oscillations  (not  resolved  in  the  figure)  are  distinctly  non- 
sinusoidal  immediately  after  the  capture.  The  amplitude  of  these 
oscillations  slowly  declines  with  time.  The  tidal  bulge  still  runs 
across  the  circumference  of  the  planet  with  a  period  equal  to  the 
orbital  period;  therefore,  the  tidal  dissipation  of  energy  goes  on. 
A  slow  shrinkage  of  the  orbit  is  the  main  source  of  tidal  energy 
for  a  planet  captured  in  spin-orbit  resonance. 

The  result  is  drastically  different  if  we  set  the  initial  conditions 
to  0(0)  =  2n,  0(0)  =  Ttl2  (Figure  5).  After  a  couple  of 
upward  swings,  the  spin  rate  jumps  through  the  resonance  point 
in  ~  80  years.  The  spin  rate  resets  abruptly  at  considerably 
lower  mean  values.  The  peak  rates  never  quite  reach  the 
resonance  value,  gradually  diminishing  with  time.  A  sidereal 
angle  0  =  ±n {2  at  M  =  0  implies  that  the  planet  is  positioned 
sidewise  with  respect  to  the  star,  that  is,  its  longer  dimension 
is  perpendicular  to  the  Sun-planet  line.  This  is  the  preferred 
configuration,  which  is  necessary  for  the  planet  to  traverse 
resonance.  The  range  of  suitable  phase  angles  at  0(0)  —  2n  and 
M(0)  =  0,  which  make  it  possible  for  the  planet  to  traverse  the 
resonance,  is  quite  small,  approximately  between  it  12  —  0.030 
and  It  12  +  0.021.  Remarkably,  at  any  other  0(0)  in  [0,  it],  the 
planet  remains  entrapped  in  this  resonance.  Thus,  the  passage 
through  resonance  cannot  occur  unless  the  planet  is  turned 
almost  sidewise  with  respect  to  the  star  at  one  of  the  periastrons. 
A  stronger  tidal  dissipation  during  the  last  circulation  before 
resonance  would  not  allow  the  planet  to  reach  this  range,  and 
capture  becomes  inevitable. 

Such  strong  requirements  for  the  orientation  of  the  planet 
at  the  point  of  resonance  may  seem  puzzling  at  first  glance. 
Flowever,  simple  qualitative  considerations  of  the  interaction  of 
the  tidal  and  triaxial  torques  explain  the  matter.  Let  us  consider 
a  triaxial  planet  in  a  resonance  20  =  (2  H-g)«,  so  that  for  each  of 
the  two  orbital  revolutions,  it  makes  exactly  2+q  rotations  about 
its  axis.  Assume  0(0)  =  0  at  an  initial  periastron  passage,  so  that 
the  longest  dimension  is  aligned  with  the  instantaneous  direction 


to  the  star.  After  two  complete  orbital  revolutions,  the  planet 
arrives  at  periastron  with  a  slight  lag  in  phase  angle  because 
the  tidal  torque  decelerated  the  planet’s  rotation.  The  planet 
therefore  is  tilted  opposite  from  the  spin  direction.  The  triaxial 
torque  will  act  to  re-align  the  planet  again,  so  that  the  average 
action  from  the  triaxial  torque  is  to  spin  up  the  planet  again. 
This  counteraction  of  the  triaxial  torque  is  symmetric,  in  that 
if  the  planet’s  spin  accelerated  during  the  two  complete  orbits 
the  torque  will  rectify  its  rotation.  Thus,  any  deviation  of  spin 
rate  from  the  resonance  value  will  be  automatically  corrected  by 
the  triaxial  torque.  The  planet  is  trapped  in  a  stable  equilibrium. 
This  is  not  the  case  when  the  initial  angle  at  periastron  is  ±7r/2. 
Any  net  lag  in  phase  angle  will  cause  a  nonzero  triaxial  torque  in 
the  same  direction,  causing  the  planet  to  continue  to  spin  down. 
In  this  configuration  of  unstable  equilibrium,  the  torque  assists 
the  tides  to  swiftly  turn  the  planet  around  and  lunge  through  the 
resonance. 

The  same  mechanics  work  for  a  circular  orbit  (e  =  0).  One 
may  ask  then,  how  a  planet  initially  at  0(0)  =  0  can  be  entrapped 
at  all,  if  at  some  point  afterward  it  inevitably  turns  sidewise 
with  respect  to  the  star?  The  answer  is  that  while  the  planet 
turns  through  it  12  with  respect  to  the  star,  the  spin  rate  will 
change  because  of  the  triaxial  torque,  and  will  no  longer  be 
equal  to  (1  +  q/2)n.  In  other  words,  the  condition  of  resonance 
passage  is  a  certain  area  of  the  two-dimensional  phase  space.^ 
Returning  to  Figure  5,  we  note  that  the  point  of  the  phase  space 
(0(0)  =  71/2,0(0)  =  2«}  does  not  actually  belong  to  this 
area,  because  the  passage  through  resonance  is  not  immediate. 
Indeed,  the  planet  has  to  realign  its  phase  space  parameters  in 
two  upward  swings  before  it  traverses  the  resonance  at  about  80 
yr  after  the  start  of  integration.  In  order  to  map  the  phase  space  of 
immediate  resonance  passages,  we  performed  several  hundred 
short-term  integrations  (for  1000  yr)  varying  the  initial  0(0)  with 
a  step  of  5  X  10“^«  and  0(0)  with  a  step  of  1  x  10“^  rad.  The 
mapped  section  of  this  space  for  our  Mercury-like  planet  and 
2;  1  resonance  is  shown  in  Figure  6  as  a  green-shaded  area.  If  the 
planet  at  periastron  turns  up  within  the  green  area,  it  immediately 
traverses  resonance.  This  “green  corridor”  is  extremely  narrow 
in  0  at  0/«  >  2  and  occupies  a  tiny  fraction  of  the  phase  space. 
The  black  dots  display  the  periastron  positions  of  the  planet. 
The  dots  are  lined  up  in  five  consecutive  libration  trajectories, 

^  In  fact,  the  relevant  phase  space  is  three-dimensional  if  we  include  the 
initial  mean  anomaly  M(0),  but  for  simplicity  of  analysis,  we  restricted  our 
study  to  the  plane  {9,  9}  by  considering  only  the  instances  of  periastron. 
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Figure  6.  Small  section  of  the  [0,6]  phase  space  around  the  2:1  spin-orbit 
resonance.  The  green-shaded  area  is  the  subspace  of  initial  parameters  from 
which  a  direct  passage  through  resonance  takes  place.  The  dots  show  consecutive 
periastron  positions  of  the  test  planet.  Parts  of  four  libration  swings  gradually 
approaching  the  point  of  resonance  are  visible  before  the  final  trajectory  hits 
the  passage  area  and  follows  it  through  resonance  into  the  domain  of  lower  spin 
rates,  0  <2n. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


each  following  one  probing  lower  spin  rates  around  0  —  it  12. 
The  final  lap  hits  the  “needle’s  eye”  and  falls  through  resonance, 
never  to  return  to  it.  Note  that  the  next  to  last  trajectory  actually 
reached  0  =  2  «  at  0  very  close  to  n  jl,  but  the  tipping  point  is 
at  a  slightly  smaller  0. 
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5.  DISCUSSION 

The  area  of  the  phase  space  above  the  resonant  spin  rate, 
through  which  a  Mercury-like  planet  can  traverse  a  resonance, 
is  very  narrow.  For  example,  the  appropriate  phase  angle  at 
periastron  for  0  —  2n  should  be  between  n 12  —  0.03065  and 
Tt 12  —  0.02575.  It  may  appear  improbable  that  the  planet,  driven 
by  the  relatively  large  triaxial  torques,  and  wandering  through  a 
good  part  of  the  phase  space  in  its  pre-resonance  evolution,  can 
hit  this  very  small  passage  opening.  And  yet,  as  we  established 
through  many  numerical  integrations,  the  test  planet  in  most 
cases  passes  resonances  higher  than  3:2,  if  the  integration  starts 
well  above  the  resonant  spin  rate.  The  camel  rarely  fails  to  go 
through  this  needle’s  eye. 

In  the  domain  of  above-resonance  spin  rates,  the  passage 
area  extends  almost  along  a  straight  line  toward  higher  0/« 
and  smaller  0.  Figure  6  shows  only  a  segment  of  this  area.  The 
needle’s  eye  area  looks  more  like  a  needle,  stretching  out  to  at 
least  0  —  Tt {2  —  0.9031,  0  jn  —  2.009,  tapering  off  to  a  point. 
On  the  other  hand,  the  lower  extents  of  libration  swings  tend 
to  probe  the  area  of  unstable  equilibrium,  i.e.,  the  minimum 
point  of  each  trajectory  is  close  to  0  =  it  12.  This  important 
fact  follows  from  the  integral  of  energy  for  free  librations  in 
the  vicinity  of  a  resonance  q  averaged  over  one  orbital  period 
(Murray  &  Dermott  2000),  which  can  be  written  as 

-  A)G20q{e)  cos(/)  =  £0,  (12) 

where  y  =  20  —  (2  +  q)M  and  Eq  is  a  constant  energy. 
When  the  rate  of  rotation  is  still  faster  than  the  resonant  value, 
20  >  {2  +  q)n,  but  it  is  at  the  minimum  of  a  libration  swing, 
0  —  min,  the  cosine  in  this  equation  should  be  equal  to  —  1, 
hence,  y  =  jr.  At  the  times  of  periastron,  the  corresponding 
value  of  0  is  7r/2.  Figures  7(a)  and  (b)  show  larger  parts  of  the 
phase  space  in  the  vicinity  of  a  2:1  resonance  passage  and  a 


Figure  7.  Phase  space  position  at  times  of  periastron  of  a  planet  passing  through 
2:1  resonance  (top)  and  captured  in  3:2  resonance  (bottom). 

3:2  resonance  capture,  respectively.  In  the  case  of  passage,  the 
trajectory  lunges  through  the  opening  depicted  in  Figure  6  from 
the  upper  zone  of  circulation  to  the  lower  zone  of  circulation, 
reversing  its  direction.  In  the  case  of  capture,  the  trajectory 
breaches  the  separatrix  between  the  upper  circulation  zone  and 
the  central  zones  of  pure  libration  where  it  follows  gradually 
tightening  loops. 

The  analogy  to  a  rotating  pendulum  by  Goldreich  &  Peale 
(1968)  can  help  us  to  qualitatively  understand  the  alignment 
required  for  the  planet  to  traverse  resonance.  The  pendulum 
is  swinging  over  the  top  of  its  support,  while  a  secular  torque 
is  slowing  down  its  rotation.  Inevitably,  the  pendulum  loses 
enough  rotational  energy  to  be  unable  to  pass  the  highest  point 
and  for  a  while  stops  close  to  the  highest  point,  after  which  it 
starts  to  rotate  in  the  opposite  direction,  but  the  secular  torque 
acts  in  the  same  direction,  this  time  assisting  the  pendulum 
in  passing  the  top  in  the  counterdirection.  Once  the  average 
y  changed  sign,  the  pendulum  continues  to  rotate  in  the 
counterdirection  with  acceleration.  This  also  explains  why  the 
trajectories  in  Figure  6  are  so  flat  at  the  minima  of  the  libration 
oscillations.  Descending  step  by  step  to  lower  spin  rates  at  each 
libration  swing,  the  trajectories  can  hardly  avoid  catching  on  the 
nearly  linear  area  of  immediate  passage  through  resonance. 

This  remarkable  fact  can  be  further  elucidated  by  analysis  of 
capture  probabilities  following  the  lines  in  Goldreich  &  Peale 
(1966).  First,  we  note  that  the  secular  tidal  torque  in  the  vicinity 
of  a  resonance  q'  (Equation  (10))  can  be  split  into  two  parts, 
one  including  the  q  —  q'  term  and  the  other  the  rest  of  the  sum. 
The  q  —  q'  term  is  an  odd  function  of  (jL>220q'  around  zero  tidal 
frequency.  The  other  term,  which  we  call  bias,  is  to  a  very  good 
approximation  constant  with  a>220q’,  because  it  is  the  sum  of  all 
other  secular  torques  at  resonances  outside  q' .  These  resonances 
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are  spaced  by  «  in  y  =  —Mrmq',  whereas  the  amplitude  of 
librations  close  to  a  resonance  is  much  smaller  (Figures  4 
and  5).  Figures  1(a)  and  (b)  show  that  the  bias  is  negative  for  the 
given  eccentricity,  implying  a  (nearly)  frequency-independent 
dissipation  of  rotation  energy.  It  is  sufficient  to  consider  two 
librations  around  the  point  of  resonance  y  =  0,  i.e.,  the  last 
libration  with  positive  y  and  the  first  libration  with  negative  y . 
Goldreich  &  Peale  (1966)  noted  that  if  the  energy  offset  from 
zero  at  the  beginning  of  the  last  libration  above  the  resonance 
is  uniformly  distributed  between  0  and  A£  =  / {T)ydt,  where 
(7)  is  the  secular  torque,  the  probability  of  capture  is 


Fcapt  — 


&E 

A£’ 


(13) 


with  5  £  being  the  total  change  of  kinetic  energy  at  the  end  of  the 
libration  below  the  resonance.  Thus,  {7)y  should  be  integrated 
over  one  cycle  of  libration  to  obtain  AZs,  and  over  two  librations 
symmetric  around  the  resonance  W220q'  =  0  to  obtain  &E.  As  & 
result,  the  odd  part  of  the  tidal  torque  at  ^  doubles  in  the 
integration  for  &E,  whereas  the  bias  vanishes.  For  the  secular 
torque  in  Equation  (10),  using  the  singular  separatrix  solution 
of  zero  energy^ 


y  —2n 


3(3  -  A) 
C 


G20q'ie) 


we  obtain 


Pcapt  — 


l  +  lTzV/  f^^W(y)dy 


cos—,  (14) 

2 


(15) 


where 


G20,(e)^.(2,  \q  -  q'\n)  Sign(q  -  q') 


W(y)  =  -  Glog:(e)K,(2,  y). 


(16) 


We  included  the  term  q  =  q'  for  bias  V  pro  forma,  because 
Kc(2,  0)  =  0.  The  integral  in  Equation  (15)  can  be  computed 
numerically  using  Equation  (14).  The  results  of  this  semi- 
analytical  estimation  of  capture  probability  as  a  function  of 
eccentricity  are  presented  in  the  top  panel  of  Figure  8.  It  is 
gratifying  to  see  that  they  are  in  agreement  with  the  results  of 
brute-force  simulations  discussed  in  Section  3.1.  At  e  =  0.205, 
the  probability  of  capture  in  a  3:2  resonance  is  1,  and  the  prob¬ 
ability  of  capture  in  a  2:1  resonance  is  approximately  0.3. 

This  fast  way  of  computing  capture  probabilities  can  be 
applied  to  exoplanets  of  terrestrial  composition,  keeping  in 
mind  that  the  probabilities  depend  on  the  degree  of  triaxiality 
through  the  parameter  (B  —  A)/ C  in  the  equation  of  separatrix 
(Equation  (14)).  Even  the  smallest  exoplanets  discovered  to  date 
tend  to  be  larger  than  the  Earth  because  of  an  observational 
selection  effect.  Larger  rocky  planets  are  likely  to  be  more 
axially  symmetric.  For  smaller  (3  —  A)/C,  keeping  all  other 
parameters  the  same,  the  curves  of  capture  probabilities  in 
Figure  8  become  steeper,  so  that  guaranteed  capture  is  achieved 
at  smaller  minimal  eccentricities.  This  is  to  be  expected  because 
if  the  triaxial  torque  is  turned  off  in  Equation  (1),  capture  is 
inevitable  at  any  resonance  where  the  tidal  torque  changes  sign. 
For  example,  the  minimal  eccentricities  of  inevitable  capture  of 
a  Mercury-like  planet  with  (3  —  A)/C  =  1.2  x  10“^  drop  to 
0.08  for  3:2,  0.20  for  2:1,  and  0.30  for  5:2  resonances. 

^  Note  that  our  y  is  twice  the  y  in  Goldreich  &  Peale  (1966). 


Figure  8.  Probability  of  capture  of  a  Mercury-like  terrestrial  planet  in  3:2,  2:1, 
and  5:2  spin-orbit  resonances.  Top:  prograde  evolution  when  the  tidal  torque  is 
generally  negative  and  the  planet  is  spun  down.  Bottom:  retrograde  evolution 
when  the  tidal  torque  is  generally  positive  and  the  planet  is  spun  up. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


Wieczorek  et  al.  (201 1)  provided  observational  evidence  that 
Mercury  was  previously  captured  in  a  synchronous  rotation.  The 
authors  of  this  paper  explain  that  this  possibility  is  achievable 
if  the  planet  started  its  evolution  with  a  retrograde  rotation.  The 
secular  tidal  torque  in  this  case  spins  the  planet  up,  first  making 
the  rotation  prograde,  and  then  driving  the  planet  toward  1:1 
resonance.  Once  captured  in  synchronous  resonance,  Mercury 
could  have  remained  there  for  an  extended  period  of  time,  until 
a  fortuitous  large  impact  drove  it  out  of  the  resonance,  abruptly 
increasing  its  rate  of  rotation.  Under  these  circumstances,  the 
planet  then  could  have  crossed  higher  resonances  in  an  upward 
order.  Since  the  previous  consideration  of  capture  probabili¬ 
ties  via  energy  balance  is  completely  symmetric  with  respect  to 
the  direction  of  y.  Figure  8,  top  panel,  is  valid  for  the  reverse 
crossing  of  resonances  for  eccentricities  below  the  equilibrium 
values.  However,  more  stringent  conditions  of  such  retrograde 
evolution  come  from  the  consideration  of  equilibrium  torques 
in  Section  3.2.  As  follows  from  Figure  3,  a  strong  impact  would 
not  be  sufficient  for  Mercury  to  leave  1 : 1  resonance.  The  orbital 
eccentricity  at  the  time  should  be  greater  than  0.29  to  over¬ 
come  the  first  barrier  of  equilibrium  torque.  Furthermore,  the 
eccentricity  should  remain  above  0.20  for  a  continuous  spin-up 
to  the  point  of  3:2  resonance.  There,  Mercury  is  guaranteed  to 
be  captured  at  any  eccentricity  between  0.20  and  0.41,  because 
the  planet  would  not  be  able  to  spin  up  any  further,  whereas 
the  probability  of  capture  at  e  >  0.41  begins  to  decline  with 
growing  eccentricity.  Figure  8,  bottom  panel,  shows  the  proba¬ 
bilities  of  capture  in  3:2,  2: 1  and  5:2  resonances  for  a  retrograde 
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evolution  of  spin  rate,  and  for  eccentricities  exceeding  the  up¬ 
per  limits  of  the  equilibrium  torque.^  Thus,  Efroimsky’s  tidal 
model  described  in  this  paper  does  not  rule  out  the  hypothesis 
by  Wieczorek  et  al.  (201 1),  but  requires,  beside  the  external  ac¬ 
tion,  fairly  high  values  of  orbital  eccentricity  during  the  ascent 
to  the  current  3:2  resonance. 

The  computations  presented  in  this  paper  for  Efroimsky’s 
model  of  tidal  torque  serendipitously  resolved  the  long-standing 
conundrum  of  Mercury’s  capture  into  3:2  resonance.  The  pre¬ 
vious  approximations  of  tidal  torque  (Goldreich  &  Peale  1968; 
Murray  &  Dermott  2000)  either  predicted  a  low  probability  of 
this  capture  with  the  current  value  of  eccentricity,  or  were  not 
able  to  reproduce  the  libration  damping.  Our  initial  computa¬ 
tions  imply  that  the  probability  of  entrapment  of  Mercury  in  3:2 
resonance  with  the  current  value  of  eccentricity  is  100%.  No 
other  outcome  would  have  been  possible  unless  the  eccentricity 
acquired  much  smaller  values  in  the  past,  or  Mercury’s  initial 
rotation  was  retrograde.  A  higher  eccentricity  likely  could  have 
resulted  in  entrapment  in  2: 1  resonance  far  in  the  past  when  the 
planet  was  approaching  this  rotation  rate. 

I  thank  the  USNO  Editorial  Board  for  helpful  suggestions 
and  a  critical  reading  of  the  original  version  of  the  paper.  Dr. 


M.  Efroimsky  inspired  this  paper  and  generously  shared  his 
insight  in  the  mechanics  of  tides. 
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Computation  of  capture  probabilities  for  retrograde  evolution  by 
Equation  (15)  is  technically  difficult  because  the  Kaula  functions  vary  rapidly 
at  high  eccentricities.  It  is  necessary  to  include  more  terms  in  the  sum  in 
Equation  (10),  e.g.,  q  =  —2,  —  1, . . . ,  8. 
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